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ABSTRACT 


Integral  equation  methods  are  described  for  calculating 
the  entire  sound  pressure  field  when  either  the  distribution  of 
velocity  or  sound  pressure  is  specified  on  an  arbitrary  closed 
surface.  The  theory  is  based  on  determining  an  equivalent  sur¬ 
face  layer  of  either  monopoles  only,  dipoles  only,  or  both  mono¬ 
poles  and  dipoles.  Appropriate  integral  equations  are  derived 
for  the  unknown  surface  monopole  and/or  dipole  density  for  each 
case  and  each  boundary  condition.  Every  closed  surface  has  twc 
infinite  sequences  of  characteristic  wave  numbers  at  each  of 
which  there  exist  an  associated  characteristic  internal  standing 
wave  and  an  associated  characteristic  external  traveling  wave 
which  satisfy  the  homogeneous  parts  of  these  integral  equations 
at  one  or  the  other  of  the  two  series  of  wave  numbers.  At 
these  wave  numbers,  and  for  particular  boundary  conditions  which 
are  specifically  derived,  all  the  integral  equations  may  have 
infinite  or  indeterminate  solutions.  The  problems  of  sound  radia¬ 
tion  by  a  pulsating  sphere  is  used  to  illustrate  the  solutions  of 
all  the  different  integral  equations  and  to  demonstrate  the  compli¬ 
cations  that  occur  at  the  characteristic  wave  numbers.  Special 
and  simple  techniques  are  described  for  approximating  each  cf 
the  integral  equations  by  a  linear  matrix  equation  with  finite 
elements  and  for  the  numerical  solution  of  the  matrix  equation. 
Special  methods  are  described  to  eliminate  the  indeterminacy  in 
the  solution  to  the  matrix  equation  near  the  characteristic  wave 
numbers. 
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INTRODUCTION 


Some  integral  equation  methods  that  are  currently  being  used  in  the 
numerical  solution  of  radiation  and  scattering  problems  are  described  and 
analyzed.  The  intention  is  not  to  provide  a  comprehensive  review  of  the 
very  extensive  literature  but  rather  to  give  a  unified  derivation  of  the 
methods  which  stresses  their  common  features  and  compares  their  relative 
merits  as  they  appear  to  the  author. 

These  methods  have  been  used  by  the  author  in  a  variety  of  specific 
engineering  problems,  e.g„,  the  distribution  of  sound  radiation  from  a 
vibrating  ship  hull,  the  directivity  pattern  of  an  array  of  sound  trans¬ 
ducers,  the  scattering  of  explosive  sound  by  the  resonant  vibrations  of 
a  ship  hull,  and  the  thrust-deduction  factor  of  a  marine  propeller.  All 
of  these  problems  required  a  solution  to  the  wave  equation  for  the  region 
outside  an  arbitrary  closed  surface,  given  some  specified  relation  between 
the  sound  pressure  and  its  normal  gradient  at  the  surface. 

However,  we  will  consider  here  only  a  restricted  version  of  this 
problem,  namely,  where  the  pressure  is  harmonic  in  time  and  where  either 
the  pressure  or  its  normal  gradient  is  specified  at  the  surface.  We  take 
the  sound  pressure  to  be  the  real  part  of  p  exp(-ikct)  and  the  gradient 
to  be  the  real  part  of  q  exp(-ikct),  where  p  and  q  are  complex.  We 
denote  an  arbitrary  point  inside,  on,  and  outside  the  closed  surface  by 
the  vectors  x,  v,  and  z.  respectively.  Hereafter  we  will  omit  the  time 
factor  from  the  equations  and  also  we  will  tacitly  assume  chat  the  real 
parts  are  taken.  Then  the  differential  equaticn  is  the  Helmholtz  equation 

(V3  +  ’3)  p(z)  *  0  (1» 

If  the  normal  gradient  of  the  pressure  is  specified  at  the  surface 
(Neumann  boundary  conditions),  the  boundary  condition  is 


dp(z) 

(3n 


irl 


q(z) 


(2) 
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where  n  is  a  unit  normal,  pointing  outward  from  the  surface  point  £.  Or 
alternatively,  if  the  surface  values  of  the  sound  pressure  are  specified 
(Dirichlet  boundary  conditions),  the  boundary  condition  is 

P(z)  “  P(£)  5  £  -  2  •  (3) 

We  also  assume  some  radiation  condition  at  infinity,  most  simply  in  the 
form 

ikz 

p(z)  «  6  z  ;  z  -►  ®  (4) 


where  z  ==  |_z| . 

Of  course,  these  simple  boundary-value  problems  occur  in  many  fields 
of  applied  physics  but  we  shall  continue  to  interpret  the  equations  solely 
in  terms  of  sound  pressure  waves.  Even  in  these  terms,  the  equations  have 
many  useful  alternative  physical  interpretations  in  sound  radiation  and 
scattering.  For  example,  there  are  two  direct  interpretations  as  radiation 
problems  (i)  The  surface  £  is  the  outer  surface  of  some  solid  body  which 
vibrates  in  an  arbitrary  patte  n  at  steady  circular  frequency  a)  =  kc, 

*  q(2)/(ipkc)  is  the  distribution  of  vibration  velocity  normal  to  the 
surface,  and  p(js)  is  the  sound  pressure  distribution  in  the  external  space; 
(ii)  the  surface  ^  is  any  geometric  surface  which  encloses  all  the  sources 
of  sound  in  the  fluid,  the  sound  sources  all  have  circular  frequency  o  =  kc, 
p(z)  is  the  external  sound  pressure  at  z,  and  q(£)/(ipkc)  is  the  normal 
component  of  the  particle  velocity  of  the  fluid  at 

There  are  also  two  straightforward  interpretations  as  simple  scatter¬ 
ing  problems  (i)  The  surface  ^  is  a  rigid  impenetrable  surface  that  scatters 
an  incoming  sound  wave  of  frequency  co  *  kc,  p(z:)  is  the  sound  pressure  of 
only  the  scattered  or  reflected  wave  (omitting  the  incident  wave),  and 
”^(z)/(lpkc)  is  the  normal  component  of  particle  velocity  in  the  incoming 
wave  at  provided  the  rigid  surface  and  its  interior  are  replaced  by 
fluid;  (ii)  The  surface  £  is  a  free  (pressure-release)  surface  that  scatters 
an  incoming  pressure  wave,  p(z)  is  the  sound  pressure  of  only  the  scattered 
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or  reflected  wave,  and  -p(^)  is  the  sound  pressure  at  ^  in  the  incoming 
wave,  provided  the  free  surface  and  its  interior  are  replaced  by  fluid. 

Finally  there  are  two  straightforward  interpretations  in  terms  of 
special  source  distributions  in  a  bounded  region  (i)  The  surface  j  is  a 
rigid  impenetrable  surface  with  a  layer  of  simple  sources.  q(”)  per  unit 
area,  placed  an  infinitesimal  distance  outside  the  surface,  and  p(£)  is 
the  external  sound  pressure;  (ii)  The  surface  £  is  a  free  (pressure- 
release)  surface  with  a  layer  of  dipoles,  n  p(v)  per  unit  area,  placed  an 
infinitesimal  distance  outside  the  surface,  and  p(£)  is  the  external 
sound  pressure.  The  last  two  interpretations  are  based  on  the  fact  (which 
we  prove  subsequently)  that  the  normal  gradient  of  the  field  of  a  layer  of 
simple  sources  and  the  field  of  a  layer  of  dipole  sources  are  each  discon¬ 
tinuous  as  the  field  point  crosses  the  layer. 

The  mathematical  methods  to  be  described  for  solving  the  boundary 
value  problems  use  integral  equations  to  determine  "equivalent1*  monopole 
and/or  dipole  layers  whose  external  field  is  the  same  as  in  the  original 
boundary  value  problem.  Hence,  we  first  review  the  properties  of  the 
fields  of  such  monopole  and  dipole  layers  and  derive  an  indefinite  number 
cf  appropriate  integral  equations.  We  thei.  show  that  each  integral  equa¬ 
tion  has  an  infinite  sequence  of  characteristic  wave  numbers,  at  which 
the  solution  becomes  infinite  or  indeterminate,  and  discuss  how  the  ori¬ 
ginal  boundary  value  problem  can  be  solved  in  such  cases.  We  then  describe 
specific  techniques  for  approximating  each  of  these  integral  equations  by 
finite  order  matrix  equations  and  for  obtaining  numerical  solutions  to 
each  of  the  matrix  equations.  These  techniques  are  specifically  designed 
for  problems  having  the  arbitrary  geometry  and  boundary  conditions  that 
commonly  occur  in  real  engineering  situations.  We  also  compare  some  of 
the  relative  merits  of  the  various  methods. 

EQUIVALENT  SURFACE  LAYERS 

The  methods  of  solution  are  based  on  the  possibility  of  finding 
some  surface  distribution  of  simple  sources  (monopoles) ,  and/or  some 
surface  distribution  of  double  sources  (dipoles),  which,  if  placed  at 
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the  position  of  the  vibrating  surface  in  an  unbounded  fluid,  would  have 
exactly  the  same  pressure  field  outside  the  surface  as  does  the  vibrating 
surface.  It  turns  out  that  there  are  an  infinite  number  of  such  equiva¬ 
lent  surface  distributions  of  monopoles  and  dipoles;*  we  shall  discuss 
only  three  special  distributions. 

The  first  equivalent  surface  distribution  is  one  having  monopoles 

only.  We  will  take  the  field  at  z  due  to  monopole  of  unit  strength  at.  y 
ikR 

to  be  e  /4jtR,  where  R  ■  £  -  y.  Then  the  field  due  to  the  entire  surface 
distribution  is 


P(z) 


fr  °(y)elkl*  •  *1 

JJ  4tt|z  -  y| 


dS 


(5) 


where  a(y)  is  the  source  strength  per  unit  area;  o(y)  must  be  specified 
so  as  to  satisfy  the  boundary  condition  of  either  Equation  (2)  or  (3), 

It  can  be  shown*  that  0(y)  must  equal  Q(y)  -  q(y),  where  q(y)  is  the 
normal  gradient  at  y  of  the  external  field  p(£),  and  Q(y)  is  the  normal 
gradient  at  y  of  an  internal  pressure  wave  P(x)  whose  surface  values 
match  those  of  the  external  field  p<£) ,  i.e.,  P(y)  =  p(y).  However, 
this  interpretation  does  not  define  3(y)  explicitly  because  the  boundary 
condition  does  not  define  both  q(y)  and  Q(y)  explicitly. 

The  second  equivalent  surface  distribution  is  made  up  only  of 
dipoles  which  are  oriented  in  a  direction  normal  to  the  surface.  Hence, 
since  we  take  the  field  at  z  due  to  unit  dipole  in  the  n  direction  at  y 
to  be  (d/3n)(e  /4«R);  where  R  ■  |£  -  >  the  field  due  to  the  entire 

surface  distribution  is 


P(*) 


eikl-  "  l\ 

_4njz  -  y| 


dS 


(6) 


1 
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where  |i(£)  is  the  local  dipole  density  per  unit  area  and  must  be  specified 
so  as  to  satisfy  the  boundary  condition  of  either  Equation  (2)  or  (3).  It 
can  be  shown^  that  p(^)  must  equal  p(^)  -  P7^).  where  p(^)  is  the  surface 
value  of  the  external  pressure  wave  p(_z)  ,  and  P;(^)  is  the  surface  value 
of  an  internal  pressure  wave  P  ;(x)  whose  normal  gradient  at  £, 
matches  that  of  p(£),  i.e.,  Q#(^)  =  q(£).  But  again  this  interpretation 
does  not  define  p(y)  explicitly. 

The  third  equivalent  surface  distribution  is  made  up  of  a  particular 
monopole  layer,  whose  surface  density  is  equal  to  the  negative  gradient 
-q(^),  and  a  particular  dipole  layer,  whose  surface  density  is  equal  to 
the  surface  pressure  p(£).  The  statement  that  the  combined  external  field 
of  these  two  particular  layers  is  equivalent  to  p(z)  is  called  the  Helmholtz 
equation 


elk! l  - 1\  n* 

ds  +  If p(*> 


a 

elk!z  "  l\ 

<§n 

-  l\m 

dS 


(7) 


1  2 

and  is  commonly  derived  from  Green's  theorem.  *  Ordinarily,  either  p(^) 
or  <l(£)  is  specified  explicitly  by  the  boundary  condition  and  the  other 
function  remains  to  be  determined. 

We  shall  call  and  n(^).  Green's  equivalent  monopole  and  dipole 

layer,  respectively,  since  they  were  apparently  first  used  by  George  Green 
in  electrostatic  problems  about  1830,  We  may  call  q(^)  and  p(y),  the 
Helmholtz  monopole  and  dipole  layers,  respectively.  If  either  a(^)  or 
n(i)  were  known,  or  if  both  p(^;)  and  q(^)  were  known,  then  p(£)  could  be 
calculated  by  simple  quadratures.  But  in  all  the  present  problems,  the 
unknown  surface  layer,  either  O(^) ,  n(^),  p(^),  or  q(^),  as  the  case  may 
be,  must  be  determined  from  the  boundary  condition  in  some  way.  In  every 
case  this  may  be  done  by  solution  of  an  appropriate  integral  equation. 

We  will  discusc  at  least  five  different  integral  equations  which  can  be 
used  for  this  purpose  when  Dirichlet  boundary  conditions  are  prescribed 
and  at  least  five  integral  equations  which  can  be  used  when  Neumann 
boundary  conditions  are  prescribed. 
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INTEGRAL  EQUATIONS  FOR  MONOPOLE  AND  DIPOLE  LAYERS 


In  mathematical  terms,  the  integral  equations  car.  all  be  interpreted 
as  limiting  forms  of  Equations  (5)  and  (6),  or  their  derivatives,  as  the 
field  point  z  approaches  a  surface  point  .  However,  the  mathematical 
processes  must  be  carefully  defined,  ar.d  restrictions  must:  be  placed  :n 
the  continuity  of  the  surface  layers  because  the  integrands,  the  integrals, 
and  the  derivatives  which  occur  are  sometimes  unbounded  and/or  discontin¬ 
uous,  A  classical,  rigorous  and  detailed  analysis  of  these  processes  was 

3 

originally  given  by  Kellogg,  His  analysis  was  confined  to  a  case  equiva¬ 
lent  to  k  s  0  in  Equations  (5)  through  (7).  But  the  additional  factor 
exp(ik]£  -  ^|)  is  only  a  minor  complication  because  this  factor  is  contin¬ 
uous  for  all  z.  Appendix  A  adapts  and  summarizes  some  of  this  analysis. 

The  arguments  are  not  meant  to  be  rigorous  but  mainly  to  give  operational 
and  detailed  meaning  to  the  abbreviated  notation  which  follows  and  which 
would  otherwise  be  ambiguous.  A  star  before  the  integral  sign  identifies 
terms  that  are  written  in  this  special  abbreviated  notation.  The  star 
may  imply  that  the  integral  is  "improper,"  that  a  "principal  value''  is 
intended,  that  some  other  special  limiting  process  is  intended,  cr  that 
some  other  special  definition  applies  which  can  be  determined  by  reference 
to  Appendix  A. 

The  integral  equations  generally  state  that  when  you  cross  a  mono- 
pole  surface  layer  the  pressure  p(r)  is  continuous,  whereas  the  normal 
gradient  q(^)  is  discontinuous  by  ar.  amount  equal  to  the  local  moncpole 
surface  density.  And,  when  you  cross  a  dipole  surface  layer,  the  normal 
gradient  is  continuous,  whereas  the  pressure  is  discontinuous  by  an 
amount  equal  to  the  local  dipole  surface  density. 

There  are  two  integral  equations  for  Green's  monopole  layer 
depending  on  the  boundary  condition.  When  q(^)  is  specified  apriori,  an 
appropriate  integral  equation  follows  from  the  discontinuity  in  the  normal 
gradient  of  a  monopole  field  and  can  be  given  in  the  form 
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where  n1  is  a  unit  normal  at  into  the  region  z.  The  integral  is  an 
abbreviated  notation  for  the  principal  value  integral  which  is  defined 
more  precisely  by  (AlO),  The  second  term  is  the  discontinuous  component 
to  1(2/)  and  changes  sign  when  ^  is  approached  from  region  x. 

When  p(^)  is  specified  apriori,  an  appropriate  integral  equation  for 
0(1)  follows  from  the  continuity  of  a  monopole  field  as  the  field  point 
passes  through  the  surface.  Hence 


-  z\ 

4*lx'  - 1\ 


dS  =  p(^') 


(9) 


where  the  integral  is  an  improper  integral  which  must  be  evaluated  as  in 

(A6). 

There  are  likewise  two  integral  equations  for  Green's  dipole  layer 

n(y),  depending  on  the  boundary  conditions.  When  q(^)  is  specified 

apriori,  an  appropriate  integral  equation  for  |i(^)  follows  from  the  con- 
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tinuity  of  the  normal  gradient  of  a  dipole  field  at  the  surface.  Thus, 
in  an  abbreviated  notation, 


da 

eik!  1  -  y\ 

dn  *dn 

.4 Ai  - 1\ 

(10) 


More  explicit  and  detailed  forms  for  this  integral  are  given  by  the 
right  hand  side  of  either  (A17)  or  (A18), 

When  p(^)  is  specified  apriori,  the  appropriate  integral  equation 
for  the  dipole  density  u(^)  follows  from  the  discontinuity  in  the  field 
of  a  dipole  layer  and  can  be  given  in  the  form 
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where  n  is  the  unit  normal  at  The  integral  term  is  an  abbreviated 
notation  for  the  principal  value  integral  which  is  defined  mon  pecis<*iv 
by  (A14).  The  second  term  is  the  discontinuous  component  t.  p(y*  ,ir.d 
changes  sign  when  the  point  is  approached  from  within  x 

There  are  at  least  three  independent  integta!  fiju.it  i  i»  I  w  i 
Helmholtz  layers  p(y;)  and  q(y),  and  each  cun  he  used  to  h  i  Diriih.c  t 
or  Neumann  boundary  conditions.  The  most  commonly  used  inu»ii<i.  «  ,  h.  -i 
is  the  Helmholtz  surface  pressure  equation,  which  Is  d<  lived  in  Apt  •  ui  i \  A 
in  Equation  (A22).  In  the  abbreviated  notation  tills  equation  is 


P(z') 

2 


dh 


An  alternative  integral  equation,  which  may  be  called  the  suifac-  . ic :  - 
equation  is  derived  in  Appendix  A,  Equation  (A25) ,  and  can  be  stated  in  chf 
abbreviated  notation  by 
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A  third  relation  which  can  be  used  as  an  integral  equation  for  the 

1  2 

Helmholtz  layers  is  the  Helmholtz  "interior"  equation  * 


'.‘kl£  -  Zl  ’ 

.  Is  -  Zl  . 


dS 


e‘Hs  -  zl 

l«-T 


where  x  is  a  field  point  inside  the  closed  surface.  The  last 
differs  from  all  the  other  Integral  equations  because  the  twe 
points  x  and  jr  do  not  have  a  common  domain. 


relatrrr. 
var iabi j 
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Finally,  we  note  the  possibility  of  obtaining  an  indefix^ite  num¬ 
ber  of  additional  integral  equations  for  either  p(^)  or  q(^)  by  a  combin 
ation  of  (12)  and  (13).  For  example,  if  we  multiply  (13)  by  \(^#)  and 
add  to  (12)  we  get 


(15) 


\{y/)  could  be  any  continuous  function  of  y7,  e.g.,  the  constant  k  \ 
However,  it  should  have  the  physical  dimensions  of  p/q,  i.e.,  of 
length,  in  order  to  maintain  dimensional  homogeneity. 

Equations  (8)  through  (15)  are  all  integral  equations  because 
the  unknown  function,  either  a(^) ,  p(y),  p(^),  or  q(^),  appears  under  an 
integral  sign.  When  the  unknown  function  also  appears  outside  the  inte 
gral,  as  in  (8)  and  (11),  and  possibly  (12),  (13),  and  (15),  the  equa¬ 
tion  is  then  an  integral  equation  of  the  second  kind,  for  which  the 

4  5 

theory  is  complete  and  well  documented.  ; 

The  properties  of  an  integral  equation  for  a  given  domain  of 
integration  are  determined  essentially  by  the  properties  of  the  kernel 
function.  There  are  four  kernel  functions  that  appear  in  all  of  these 
integral  equations.  They  all  have  simple  physical  interpretations,  and 
they  are  all  simply  related  to  each  other.  These  kernel  functions  are 
defined  as  follows  in  the  abbreviated  notation  as  well  as  in  a  detailed 
notation,  with  R  *  ^ ;  R  |  r|  ;  R  =  R/R. 


-  l\  eikR 


(16) 


is  the  field  at  due  to  a  unit  monopole  at  ^  in  a  free  unbounded  field 
and  is  a  function  of  the  two  vector  field  points. 
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is  the  field  at  ’  due  to  a  unit  dipole  n  at  £  in  an  unbounded  field  and 
is  a  function  of  the  two  positions  and  the  direction  n.  The  subscript  2 
indicates  that  a  derivative  is  taken  with  respect  to  the  second  variable. 
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is  the  field  gradient  in  the  n/  direction  at  due  to  a  unit  monopale 
at  ^  in  an  unbounded  field  and  is  a  function  of  the  two  points  and  the 
direction  n'. 
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(19) 


is  the  gradient  in  the  n‘  direction  of  the  field  at  yj  due  to  a  unit 
dipole  n  at  £  in  an  unbounded  fluid  and  depends  on  the  two  points  and 
the  two  directions. 

These  kernels  all  have  weak  singularities  at  ^  ^  but  3  prin¬ 

cipal  value  integral  exists  in  all  cases.  Note  that  the  two  functions 
M  and  E  are  each  self-symmetric  (but  not  complex-conjugate-symmetric) 


11 


to  an  interchange  in  the  two  field  points  and  the  two  directions,  i.e., 

*  M(v,^ ;)  and  E(x  Vfl/j.x/ji)  *  Vn).  However,  the  two  func¬ 

tions  D  and  N  are  mutually  symmetric  (but  not  ccmplex-con jugate-symmetric! 
to  an  interchange  in  the  same  variables,  i.e.,  D(^,^ ' /n  ')  =  N(^  Vn  '  f'£) . 

The  integral  equations  for  these  equivalent  monopole  and  dipole 
layers  may  have  no  finite  solution,  or  possibly  an  indeterminate  solution, 
under  certain  commonly  occurring  practical  conditions.  This  may  seem 
surprising  because  the  boundary  value  problem,  for  the  space  external  to 

an  arbitrary  surface,  posed  by  Equations  (1)  through  (4)  always  has  a 

2 

solution  and  the  solution  is  unique.  However,  the  same  integral  equa¬ 
tions  can  also  be  intermediate  steps  in  solving  a  boundary  value  problem 
for  the  region  interior  to  the  surface,  and  it  should  be  obvious  on 
physical  grounds  that  the  solution  to  the  interior  problem  becomes  inde¬ 
terminate  at  the  resonance  frequencies  of  interior  standing  waves.  Accord 
ingly,  we  will  discuss  in  the  next  section  the  possible  standing-wave 
solutions  in  the  space  interior  to  any  arbitrary  closed  surface. 

INTERNAL  STANDING  WAVES  AND  ASSOCIATED 
EXTERNAL  TRAVELING  WAVES 

We  shall  assume  without  proof  the  assertion^  that  for  every  c3.osed 

surface,  an  unbounded  sequence  of  wave  numbers  exists  at  which  standing 

waves  can  e:-ist  inside  the  surface  with  any  prescribed  combination  of 

field  and  normal  gradient,  of  (^)  4-  &Q(i)  *  0,  at  the  surface  where  C*  and 

$  are  real  constants.  We  consider,  in  particular,  the  two  independent 

series  of  characteristic  wave  numbers,  those  for  which  {5=0  and  P(^)  =  0 

and  those  for  which  OfsO  and  Q(%)  ■  0.  We  denote  the  characteristic  wave 

numbers  of  the  first  kind  by  l,;  i  =1,2,...  the  interior  standing  wave 

C1)  1 

by  (x)  (normalized  in  any  convenient  way),  and  its  normal  gradient 

at  the  surface  by  Q^^).  Also  we  define  an  associated  traveling  wave 

(  x )  * 

p^  (z)  for  the  external  space  by  the  condition  that  the  normal  gradients 
match  at  i.e. , 
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Likewise,  we  denote  the  characteristic  wave  numbers  of  the  second 


and 


kind  by  m  and  the  internal  normalized  standing  wave  by  P'  ' (x) , 

(a)  1  ~ 

define  an  associated  traveling  wave  (z)  for  the  external  space  by  the 

condition  that  the  fields  match  at  i.e., 

=  p*8)(z>  (i 


Note  that  for  the  inteiior  of  any  closed  surface,  the  fundamental 
characteristic  standing  wave  of  the  seconc  class  is  (x)  =  constant 
with  k  =  =  0i  For  this  wave  function  clearly  satisfies  the  Helmholtz 

equation  and  the  boundary  condition  '  =  0. 

More  generally,  for  any  surface  shape,  the  characteristic  wave 
number  may  be  degenerate,  i.e,,  there  may  exist  a  finite  number  of 
independent  wave  functions  P^x)  and  p ^ j C^.) ;  j  ■  l,2,...jmax,  whe*e  J 
depends  on  the  surface  shape  and  wave  number.  However,  we  shall  omit 
any  further  reference  to  the  degeneracy  of  the  wave  functions  in  order 
to  simplify  the  notation. 

We  show  in  Appendix  B  that  the  equivalent  monopole  and  dipole 
layers  for  these  four  series  of  characteristic  wave  functions  have  seme 
simple  and  remarkable  properties.  Thus  for  characteristic  functions  of 
the  first  kind,  i.e.,  if  k  *  A^,  at  which  there  can  exist  internal  stand¬ 
ing  waves  with  zero  pressure  at  the  surface,  then 
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For  characteristic  functions  of  the  second  kind,  i.e. ,  if  k  =*  m^,  (at 
which  there  can  exist  internal  standing  waves  with  zero  normal  gradient 
at  the  surface) 
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Note  that  the  integrals  in  (22a),  (22c),  (23a),  (23c),  (26a),  (26c),  (27a), 
and  (27c)  need  not  be  starred. 

Equation  (22c)  shows  that  for  an  arbitrary  closed  surface  there 

exists  one  sequence  of  wave  numbers  (the  characteristic  wave  numbers  of 

the  first  kind)  at  which  an  associated  surface  distribution  of  monopoles, 

q^(^)  produces  no  external  field,  and  (27c)  shows  that  there  is  a 

second  sequence  of  wave  numbers  at  which  an  associated  surface  distribution 
(a) 

of  normal  dipoles  p:  /(^)  produces  no  external  field.  These  results 

1  7 

generalize  a  remark  of  Lamb  who  showed  that  when  sin  ka  =  0,  "a  uniform 
distribution  of  simple  sources  over  a  sphere  of  radius  a  produces  no 
effect  at  external  points,"  Thus  the  surface  distributions  that  produce 
a  null  external  field  exist  for  every  surface  shape,  for  dipole  distribu¬ 
tions  as  well  as  simple  source  distributions,  and  they  exist  in  an  infinite 
number  of  characteristic  patterns.  In  fact  at  k  *  0,  a  uniform  distribution 
of  normal  dipoles  over  any  closed  surface  produces  no  external  field. 

However,  these  results  are  principally  of  consequence  in  the  solu¬ 
tion  of  the  integral  equations  because  (22b)  and  (29)  demonstrate  a 
family  of  solutions  to  the  homogeneous  parts  of  the  integral  equations 
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of  th'i  first  kind,  whereas  (23b),  (24),  (27b),  and  (28)  demonstrate  a 
family  of  solutions  to  the  homogeneous  parts  of  the  integral  equations  of 
the  second  kind.  And  it  is  well-known  that  the  solution  to  any  linear 
integral  equation  becomes  indeterminate  under  conditions  in  which  the 
homogeneous  part  of  the  equation  has  a  solution. 

EXISTENCE  AND  UNIQUENESS 

The  mathematical  condition  for  the  existence  of  a  unique  solution 
to  any  linear  integral  equation  of  the  second  kind  is  known  as  the  Fredholm 
Alternative.  As  applied  to  (8),  for  example,  this  states  that  either 
the  integral  equation  with  a  specified  k  and  a  specified  surface  shape  has 
a  unique  and  regular  solution  for  all.  choices  cf  the  known  function  q(^)  or 
the  homogeneous  ,/art  of  the  integral  equation  has  a  rontrivial  solution  and 
the  adjoint  homogeneous  integral  equation  has  a  nontrivial  solution.  In  the 
alternative  case  the  integral  equation  has  a  regular  (but  not  unique)  solu- 
tion  only  when  the  known  function  is  orthogonal  to  every  solution  of  the 
adjoint  homogeneous  equation. 

In  the  present  context,  two  functions  fx(^)  and  fa  (^)  are  ortho¬ 
gonal  if 

||  Ms  dS  =  ||  fits  dS  -  0  (30) 

where  the  "overbar"  denotes  the  complex  conjugate.  And  by  the  adjoint 
homogeneous  equation  we  mean  a  homogeneous  integral  equation  of  the  second 
kind  whose  kernel  function  is  obtained  from  the  original  kernel  by  inter¬ 
changing  the  two  variables  and  taking  the  complex  conjugate.  Thus  adjoint 
homogeneous  equations  to:  (i)  Equation  (8)  for  a(jr),  (ii)  (11)  for  p(£), 
(iii)  (12)  for  p(^),  and  (iv)  (13)  for  q(j;)  are  the  complex  conjugates  of 
(i)  (23b),  (ii)  (27b),  (iii)  (24),  and  (iv)  (28),  respectively. 

The  special  relations  of  (22)  through  (29)  demonstrate  that  solu¬ 
tions  to  the  homogeneous  integral  equations  for  sound  radiation  do 
occur  and  that  they  occur  at  the  characteristic  wave  numbers  of  the 
first  and  second  kind  i ^  and  m^.  Also,  they  identify  the  solutions  of 
the  homogeneous  equations  to  be  the  surface  values  of  the  characteristic 
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/^\  /g\ 

traveling  waves  of  the  first  and  second  kind--p;'  '(y),  q)  '  (%)  >  p)  '(%.)> 

(2)  111 

and  q^  '  (y) .  Furthermore,  we  show  in  Appendix  B  that  solutions  to  these 
homogeneous  integral  equations,  excluding  (14)  and  (15),  can  occur  only 
at  the  characteristic  wave  numbers  or  m^. 

Accordingly,  we  list  below  the  conditions  for  the  existence  and 
uniqueness  of  solutions  to  the  integral  equations  for  c(^),  n(y),  p(_y) , 
and  q(^).  These  are  stated  as  conditions  on  the  prescribed  wave  number  k 
and  on  the  prescribed  boundary  condition  p(^)  or  q(^).  The  conditions 
may  be  verified  by  systematic  use  of  the  Fredholm  Alternative  and  by 
appeal  to  the  special  relations  of  (22)  to  (29). 

1.  A  unique  equivalent  monopole  layer  O(y)  exists  as  a  solution 

to  either  integral  equation  (8)  or  (9)  unless  k  is  equal  to  a  character¬ 
istic  wave  number  of  the  first  kind,  If  k  =  and  if  either 

JJ  qp^dS  =  0  or  JJ  pcf^dS  -  0,  see  Equation  (B21);  then  a  particular 
solution  to  either  (8)  or  (9)  for  O(^)  exists.  Any  finite  multiple  of 
q^(2)>  however,  can  be  added  to  O(^)  with  no  effect  on  the  external 
field  as  calculated  from  (5).  If  k  *  and  JJ  qp^dS  =  Jj  ptf^dS  ^  0, 
no  finite  J(^)  exists. 

2.  A  unique  equivalent  dipole  layer  n(jjr)  exists  as  a  solution  to 

either  integral  equation  (10)  or  (11)  unless  k  is  equal  to  a  characteristic 

wave  number  of  the  second  kind,  m. .  If  k  =  m.  and  if  either  |[  pP^dS  =  0 

fT  _  (2)  *  *  JJ  i 

or,  equivalently,  if  I  qp)  dS  =  0,  then  a  particular  solution  n(y) 

1  (a) 

exists.  However,  any  multiple  of  p^  (^)  can  be  added  to  ^i(^)  with  no 
effect  on  the  external  field  as  calculated  from  (6).  If  k  =  m^  and 

JJ  ptl^dS  =  JJ  qp^dS  ^  0,  then  no  finite  |i(^)  exists. 

3.  A  particular  solution  to  the  Helmholtz  surface-pressure  equation, 
(12),  for  either  p(^r)  or  q(^),  exists  for  all  k  and  is  unique  unless  k  =  l. 

(i)  1 

When  k  *»  j^,  then  any  finite  multiple  of  q£  can  be  added  to  a  particular 

solution  for  q(^)  with  no  effect  on  the  external  field  as  calculated  from 

(6).  But  also  when  k  =  A,  the  particular  solution  to  (12)  for  p(j£)  is 

1  (l) 

indeterminate  by  an  arbitrary  multiple  of  p^  (y),  and  the  corresponding 
formal  solution  of  (6)  for  p(^)  is  indeterminate  by  the  same  multiple  of 
p^(z).  The  correct  multiple  of  p^(^)  which  fits  the  specified 'boundary- 
condition  mu3t  be  determined  by  some  supplementary  condition. 
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4.  A  particular  solution  to  the  surface-gradient  equation  (13),  for 

either  p(y)  or  q(y),  exists  for  all  k  and  is  unique  unless  k  =  m. .  When 

^  (a)  * 

k  =  m  ,  any  finite  multiple  of  p)  7  (y)  can  be  added  to  the  particular  solu¬ 
tion  of  (13)  for  p(y)  with  no  effect  on  the  external  field  p(jz)  as  calculated 
from  (6).  But  also  when  k  =  m  the  particular  solution  to  Equation  (13)  for 

^  /jg  \ 

q(y)  is  indeterminate  by  an  arbitrary  multiple  of  q)  7  (y) ,  and  the  consequent 

1  (a) 

solution  of  (6)  for  p(£)  is  indeterminate  by  the  same  multiple  of  p'  (z). 

(3}  1 

Again  the  correct  multiple  of  p^  * ( z )  which  fits  the  specified  boundary 
condition  may  be  determined  by  some  supplementary  condition. 

5.  The  existence  and  uniqueness  of  solutions  to  the  Helmholtz  interior 
equation  (14)  depend  upon  the  position  of  the  interior  points  x  as  well  as  on 
the  wave  number  k.  We  can  list  only  some  of  the  conditions  under  which  the 
solution  of  (14)  is  indeterminate.  These  occur  when  k  =  l,  or  m.,  and  all 

f  X  *  X  /  \ 

x  are  on  a  nodal  surface  of  the  internal  standing  waves  p)  7 (x)  or  ’ (x) , 

respectively.  If  k  =  l  ,  and  p)  7 (x)  =  0,  then  any  finite  multiple  of 

(i)  li- 

q  (y)  can  be  added  to  the  particular  solution  of  q(y)  with  no  effect  on 

the  external  field  as  calculated  from  (6),  Under  the  same  conditions  the 

particular  solution  to  (14)  for  p(y)  is  indeterminate  by  an  arbitrary 

multiple  of  P^(y),  and  the  consequent  solution  of  (6)  for  p(z)  is  indeter- 

minate  by  the  same  multiple  of  p^  •  (£),  If  k  ■  and  P^  7 (x)  =  0,  the 

particular  solution  to  Equation  f  14)  for  q,  (y)  is  indeterminate  by  an 

(a)  1  " 

arbitrary  multiple  of  q'  (y) ,  and  the  consequent  solution  to  (6)  for  p(z)  is 
indeterminate  by  the  same  multiple  of  pj*^(z).  Under  the  same  conditions 

(a)  1 

(k  =  m.  and  P:  (x)  =  0)  the  particular  solution  to  (14)  for  p(y)  is 

*■  (a) 

indeterminate  by  an  arbitrary  multiple  of  p^  (y);  ideally  this  indeter¬ 
minate  component  has  no  external  field. 

6.  The  existence  and  uniqueness  of  solutions  to  the  combination 
equation  (15)  depend  upon  the  choice  of  the  function  \(v)  as  well  as  on 
the  wave  number  k.  However,  if  \(y)  is  chosen  as  a  positive  constant 
independent  of  y,  it  seems  plausible  that  the  integral  equation  would 
have  an  infinite  sequence  of  characteristic  wave  numbers  with  values  that 
interlace  the  and  interlace  the  m^. 
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7,  The  conditions  for  the  existence  and  uniqueness  of  ary  equiva.tnt 
monopole  layer  to  be  used  in  Equation  (5)  are  the  same  as  the  conditions 
for  the  solutions  of  integral  equations  (8)  and  (9).  Likewise  the  condi¬ 
tions  for  an  equivalent  dipole  *ayer  to  be  used  in  (6)  are  the  same  as  the 
conditions  for  the  solutions  of  integral  equations  (10)  and  (11).  However, 
there  always  exist  correct  and  unique  Helmholtz  layers  p(y)  and  qvy)  which 
can  be  used  in  (7),  even  though  the  solutions  to  one  of  the  integral  equa¬ 
tions  (12)  through  (15)  may  be  indeterminate. 

In  principle,  the  characteristic  wave  numbers  for  any  one  of  thes? 
integral  equations  are  discrete  line  values;  i.e.,  the  singularities  in  tnt 
solutions  would  occur  only  if  k  were  exactly  equal  to  or  and  in 
the  case  of  (14),  only  if  x  were  exactly  on  a  nodal  surface.  However ,  in 
practice,  the  integral  equations  are  solved  by  numerical  methods  which 
necessarily  generate  an  indeterminate  error,  depending  on  the  precision  of. 
the  computing  machine  and  the  nature  of  the  computation  algorithms.  Thus, 
there  will  be  a  band  cf  values  for  k  surrounding  each  and  m^,  within 
which  the  solution  to  the  integral  equation  has  an  indeterminate  component 
or  a  singularity;  for  (14)  there  will  be  a  layer  of  positions  about  each 
nodal  surface.  Furthermore,  for  every  surface  shape,  the  average  spacing 

-3 

in  k  between  successive  or  m^  decreases  as  k  ,  and,  eventually,  the 
bands  (or  layers)  must  overlap.  Thus,  for  every  surface  shape  and  every 
integral  equation,  there  is  a  value  of  k  beyond  which  k  is  arbitrarily 
close  to  a  characteristic  value  and  beyond  which  a  numerical  solution  of 
an  integral  equation  will  generate  an  indeterminate  error. 

In  summary,  the  numerical  solutions  for  p(x)  or  q(^)  will  be 
indeterminate  by  an  arbitrary  multiple  of  some  particular  surface  distribu¬ 
tion  when  the  wave  number  k  is  too  large  or  when  k  is  too  close  to  some 
characteristic  value.  Under  similar  conditions  on  k,  a  finite  solution 
for  a(£)  or  ti(^)  will  exist  only  for  special  boundary  conditions,  even 
then,  the  indeterminacy  will  occur. 

However,  these  difficulties  do  not  mean  that  the  integral  equations 
cannot  be  used  in  the  neighborhood  of  these  characteristic  wave  numbers. 

In  some  cases,  e.g. ,  the  solution  of  (12)  for  q(^)  or  the  solution  of  (13) 
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for  p(^),  the  indeterminate  component  in  the  solution,  has  no  external 

field  and,  hence,  will  not  cause  an  error  in  the  subsequent  calculation  of 

the  external  field  from  either  (5),  (6),  or  (7).  In  other  cases  where  the 

8 

indeterminate  component  does  have  an  external  field,  it  has  been  possible 
to  supplement  the  original  integral  equation  with  enough  additional  condi¬ 
tions  to  fix  the  correct  amplitude  of  this  component.  Even  in  those  cases 
where  no  finite  solution  exists,  the  experience  described  in  the  next 
section  suggests  that  it  may  be  possible  to  factor  out  the  infinite  compon¬ 
ent  from  the  external  field. 

UNIFORMLY  PULSATING  SPHERE 

To  illustrate  all  of  these  methods,  and  to  demonstrate  the  complica¬ 
tions  which  occur  at  the  resonance  wave  numbers,  we  will  consider  the  very 
simple  problem  of  a  sphere  of  radius  a,  pulsating  with  uniform  radial 
velocity  in  a  uniform  unbounded  medium.  Assume  that  either  the  surface 
pressure  pQ  exp(-ioit)  or  the  normal  gradient  at  the  surface  qQ  exp(-ioit) 
is  specified  apriori.  Then  explicit  analytical  solutions  to  any  one  of 
the  integral  equations  (8)  through  (14)  can  be  obtained  because  all  the 
unknown  functions  are  constants  which  can  be  factored  outside  the  integral 
sign  and  because  all  of  the  integrals  can  be  reduced  to  functions  of  the 
single  variable  R  *  J r  -  y|  by  taking  the  surface  element  in  the  shape  of 
a  spherical  zone  between  the  two  circles  at  distances  R  and  R  +  dR  from 
the  field  point  r.  The  method  of  solution  is  demonstrated  in  Appendix  C. 

By  these  methods  the  explicit  solutions  to  (8)  and  (9)  for  the 
equivalent  monopole  density  are  readily  found  to  be 

°<z'>  *  ikT^r  %  <31> 

-Gets)  "o  <32> 

Likewise,  the  solutions  to  (10)  and  (11)  for  the  equivalent  dipole 
densities  are 
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A  formal  solution  to  the  Helmholtz  surface  equation  (12)  for  either  of 
the  Helmholtz  densities  p(^)  or  q(v)  is 

(sin  ka)  q(^r')  =  (sin  ka)  (ik  -  i)  p(^')  (35' 

A  formal  solution  to  the  surface  gradient  equation  (13),  for  either  of 
the  two  Helmholtz  densities  is 

(ka  -  tan  ka)  q(^y)  ■  (ka  -  tan  ka)  (ik  -  — )  p(y;)  36) 

andjfinally,  a  formal  solution  to  the  Helmholtz  interior  equation  (14) 
for  either  of  the  two  Helmholtz  layers  is 

(sin  kx)  q(jr#)  =  (sin  kx)  (ik  -  ^)  p(^')  (37) 

Hence  in  this  problem,  every  one  of  the  integral  equations  has  a 

formal  solution  for  all  values  of  k.  However,  each  solution  has  some 

singularity  at  those  special  values  of  k  where  the  homogeneous  part  cf 

the  integral  equation  would  have  a  solution.  Thus,  the  formal  solutions 

for  o(^)  become  infinite  where  sin  ka  =  0,  and  this  condition  identifes 

the  characteristic  wave  numbers  of  the  first  kind  l  ,  where  there  can 

i' 

exist  internal  standing  waves  with  zero  surface  pressure.  The  formal 
solutions  for  p(;jr)  become  infinite  where  tan  ka  =  ka,  and  these  identif> 
the  characteristic  wave  numbers  of  the  second  kind  where  there  can  exist 
internal  standing  waves  with  zero  surface  gradient.  The  solutions  to 
the  Helmholtz  surface  pressure  equation  (12)  become  formally  indeterminate 
when  sin  ka  =  0,  again  at  the  characteristic  wave  numbers  of  the  first 
kind.  However,  the  solutions  to  the  surface  gradient  equation  (13'  for  thr 
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s^aie  Helmholtz  layers  become  indeterminate  at  the  characteristic  wave 
numbers  of  the  second  kino.  And  the  solutions  tc.  the  Helmholtz  interior 
equation  are  formally  indeterminate  when  sin  kx  =  0.  It  is  interesting 
that  this  can  occur  not  only  at  the  characteristic  wave  numbers  l and 
m^  but  also  at  any  wave  number  k  2  rr/a. 

Tf  we  disregard  the  apparent  singularities  and  simply  substitute 
either  (31)  or  (32)  into  (5),  or  substitute  (33)  or  (34)  into  (6),  or 
substitute  either  (35) ,  (36) ,  or  (37)  into  (7),  we  will  get  expressions 
for  the  external  field  p(z)  which  are  formally  indeterminate  at  the 
critical  wave  numbers.  But  in  each  case  we  can  simply  divide  out  the 
indeterminate  factor,  or,  more  generally,  we  can  require  that  the  solu¬ 
tion  for  p(z)  be  a  continuous  function  of  k.  In  this  way  we  obtain  the 
correct:  and  unique  external  field  pressure. 

,  x  V  a  ik(z  -  a)  1 

■  te;—  ;■ 


K.  a  ik(z  •*  a) 

z)  =  p  -  e 
o  7 


which  is  valid  for  all  k,  even  at  those  critical  wave  numbers  where 
either  a  finite  equivalent  monopole  or  dipole  layer  does  not  exist. 


COMPARISON  OF  INTEGRAL  EQUATION  METHODS 

We  have  derived  at  least  five  different  integral  equation  methods 
for  solving  the  Dirichlet  problem  and  at  least  five  different  integral 
equation  methods  for  solving  the  Neumann  problem.  Equation  (10)  for 
the  Neumann  problem  and  (13)  and  (15)  for  either  the  Dirichlet  or 
Neumann  problems  appear  to  be  new  in  sound  radiation  and  scattering 
problems;  however,  all  of  these  methods  have  a  long  and  continuous 
history  of  application,  development,  and  rediscovery  in  problems  in  fluid 
flow,  elasticity,  acoustic  waves,  and  electromagnetic  waves.  Mikhlin'* 
cites  a  long  bibliography  of  older  publications  by  Soviet  authors,  and 
there  is  no  doubt  that  an  equally  long  bibliography  can  be  prepared 


with  non-Soviet  authors.  However,  the  great  power  of  these  methods  could 
not  be  realized  until  the  development  of  high  speed  computers,  and  we  par 
ticularly  wish  to  stress  the  more  recent  methods  and  applications  using  these 
computers.  Practical  experience  in  such  applications  to  acoustic  problems 

seems  to  be  confined  to  solutions  of  the  Neumann  problem,  most  often  b> 

9-13  14-19 

calculating  a(^)  from  (8)  or  by  calculating  p(^)  from  (12)-; 

however,  there  is  some  limited  experience  with  solutions  of  the  Helmho.tz 
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interior  equation  (14)  for  p(^). 

It  seems  generally  preferable  to  base  a  solution  of  either  the 
Dirichlet  or  Neumann  problems  on  a  calculation  of  the  unknown  q(^,  cr  p 
rather  than  on  a  calculation  of  the  equivalent  monopole  layer  o(y'  cr  the 
equivalent  dipole  layer  |i(2)»  This  is  because:  (i)  the  functions  p(y'  or 
q(z)  have  immediate  physical  significance  and  may  even  be  determined  exper¬ 
imentally;  (ii)  in  many  practical  applications,  e,g. ,  low-frequency  radia¬ 
tion  by  longitudinal  vibrations  of  a  slender  body,  the  Helmholtz  integral 
term  in  p(^)  contributes  negligibly  to  the  far  field,  and  there  is  no 
need  to  solve  an  integral  equation  to  determine  p(y)»  (iii)  in  cases  of 
high-frequency  radiation  by  an  arbitrary  body  to  the  far  field  there  is  a 
simple  approximation  to  the  unknown  function  p(y)  or  q(y)  which  again  makes 
it  unnecessary  to  solve  an  integral  equation;  and  (iv)  when  the  wave  number 
k  is  too  large  or  too  close  to  a  characteristic  wave  number,  a  finite  r(v) 
or  p(^)  may  not  exist  but  a  unique  and  finite  p(^)  or  q(^)  always  exists. 

However,  the  relative  merits  of  integral  equations  cf  the  first 
kind  and  of  the  second  kind  are  not  so  easy  to  assess,  Equaticns  (8' 
and  (12)  are  both  equations  of  the  second  kind  for  the  Neumann  problem, 
having  similar  kernel  functions  and  with  the  same  spectrum  of  critical  wave 
numbers.  Equations  (11)  and  (13)  are  both  integral  equations  of  the  secend 
kind  for  the  Dirichlet  problem,  again  having  similar  kernels  and  witr.  ;he 
same  spectrum  of  critical  wave  numbers.  Ordinarily,  we  would  expect  that 
the  equations  of  the  second  kind  are  preferable  and  easier  to  solve  than 
the  first  kind  because  the  theory  of  the  former  is  complete  and  well  d  xu- 
mented,  whereas  it  is  still  possible  that  unexpected  difficulties  can  arise 
in  solutions  of  the  latter  equations.  However,  the  equations  of  the  first 
kind  may  have  some  special  advantages  at  the  critical  wave  numbers.  For 
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example,  if  we  solve  the  Dirichlet  problem  at  k  =  by  using  (12),  an 
equation  of  the  first  kind,  then  the  solution  is  indeterminate  by  some 
multiple  of  However,  this  indeterminate  component  has  no  effect  on 

the  external  field.  Likewise  a  solution  to  the  Neumann  problem  at  k  =  nu  , 
obtained  by  using  an  equation  of  the  first  kind  (13),  has  an  indeterminate 

(a) 

component  p^  '  without  an  external  field. 

It  may  sometimes  be.  useful  to  seek  a  simultaneous  solution  of  two  of 
these  integral  equations  rather  than  only  one.  For  example,  one  method 
fer  solving  the  Neumann  problem  when  k  ~  is  by  a  simultaneous  solution 
of  the  Helmholtz  surface  equation  (121  and  the  Helmholtz  interior  equation 
(14).  The  latter  must  be  satisfied  at  a  sufficient  number  of  interior 
points  x  to  fix  the  amplitude  of  any  indeterminate  component  to  the  solu¬ 
tion  of  (12).  The  scheme  has  value  only  when  the  wave  number  k  is  close 
enough  to  a  characteristic  wave  number  of  the  first  kind  l  to  generate 
an  indeterminate  component  in  the  solution  to  (12)  but  when  the  interior 
points  x  can  be  picked  not  too  close  to  the  nodal  surface  of  an  internal 
standing  wave.  When  k  is  sufficiently  large,  the  layers  of  critical 
positions  about  each  nodal  surface  will  overlap,  and  this  method  cannot 
be  used. 

An  alternative,  and  preferable,  scheme  is  to  solve  either  the 
Neumann  or  Dirichlet  problems  by  a  simultaneous  solution  of  both  surface 
integral  equations  (12)  and  (13),  The  critical  wave  numt srs  for  these 
equations  are  the  characteristic  wave  numbers  of  the  first  and  second 
kind,  respectively.  They  do  not  coincide  except  for  special  surface 
shapes,  e.g. ,  a  rectangular  box,  and  then  the  characteristic  wave  functions 
are  orthogonal  on  the  surface  and  need  cause  no  difficulty.  Nevertheless, 
it  is  likely  that  any  algorithm,  and  Equation  (15)  can  be  viewed  as  one 
such  algorithm,  which  "selects"  the  "accurate"  solution  from  an  over¬ 
determined  combination  of  two  integral  equations,  will  tend  to  introduce 
its  own  set  of  characteristic  wave  numbers.  In  any  case  the  numerical 
calculation  will  undoubtedly  become  ill-conditioned  at  sufficiently  high 
wave  numbers  where  the  critical  bindwidths  about  l  and  m^  overlap. 
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REDUCTION  OF  INTEGRAL  EQUATIONS  TO  MATRIX  EQUATIONS 


The  first  step  in  obtaining  a  numerical  solution  to  any  of  these 

integral  equations  for  either  the  Neumann  or  Dirichlet  problem  is  to  con 

vert  the  equation  into  a  matrix  equation,  i.e.,  into  a  set  of  simultaneous 

algebraic  equations.  The  procedures  that  we  will  describe  are  designed  tot 

the  most  general  case,  where  the  surface  has  some  arbitrary,  nonanaiytical 

shape,  and  where  the  known  boundary  condition  in  q(^)  or  p(^)  is  likewise 

arbitrary  and  nonanaiytical.  Even  in  those  cases  where  the  surface  is  a 

spheroid  or  an  ellipsoid,  and  where  in  principle  the  kernel  function  can 

be  expanded  in  terms  of  separable  wave  functions,  it  may  be  more  accurate 

to  solve  the  integral  equation  by  these  general,  nonanaiytical,  techniques 

because  the  higher  order  wave  functions  are  not  tabulated  and  are  too 

1 7 

difficult  to  compute. 

The  basic  procedure  is  to  approximate  the  integral  equation  that  is 
specified  over  a  continuous  range  of  positions  on  the  surface  by  a  matrix 
equation  that  is  specified  at  only  a  finite  number  of  stations  on  the  sur¬ 
face.  To  do  this  we  first  pick  a  finite  number  N  of  representative  stations 
^5  j  “  1,2,..,N  on  the  surface  and  then  associate  each  station  with  an 
element  of  area  Sj,  a  normal  direction  n j ,  and  some  representative  value 
Cj  for  the  local  radius  of  curvature  of  the  surface.  The  stations  need 

not  be  spaced  with  uniform  density,  and  the  c  ements  of  surface  area  need 

18 

not  have  any  special  or  uniform  shape  The  density  of  stations  on  the 
surface  should  be  high  in  regions  of  high  surface  curvature  and  in  regions 
where  p(^)  or  q(v)  changes  most  rapidly,  except  that  no  station  should  be 
placed  at  a  discontinuity  of  p(^)  or  q(£).  The  shape  of  each  surface  ele¬ 
ment  need  not  be  specified;  the  magnitude  of  each  area  Sj  should  be  deter¬ 
mined  by  distributing  the  total  surface  area  E  Sj  among  the  N  elements  in 
approximately  inverse  proportion  to  the  local  station  density  and  in 
accordance  with  any  convenient  bookkeeping  scheme.  In  practice,  it  may 
be  convenient  to  assume  that  the  elements  have  some  particular  shape  e,g,  , 
quadrilateral  or  triangular;  however,  the  shape  has  no  effect  on  the 


calculations.  The  general  prescription  is  that  p(^),  q(2)>  the  direction 
of  the  normal  n(^),  and  the  curvature  of  the  surface  should  each  be  reason 


ably  uniform  over  a  compact  area  Sj  about  the  point 
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We  now  approximate  each  of  the  integrals  that  occurs  in  (8) 
through  (14)  by  a  sum  of  terms  evaluated  only  at  the  N  stations  For 

example,  if  the  unknown  function  is  and  the  kernel  function  is 


we  sec 


N 

JJ  Hi)  Hlyl)  d9  wjk  f(ik)  K 

k=i 


Jk 


(40) 


for  i  =  1,2,...N.  The  weighting  factor  W.,  depends  upon  the  quadrature 
formula  being  used.  In  the  simplest  and  most  common  case  W^k  =  Sk{  how¬ 
ever,  for  special  surface  shapes  or  for  special  ways  of  ordering  the 
points  2.j;  it  is  possible  to  make  a  more  sophisticated  choice  for  the 
weighting  factor.  The  value  for  the  kernel  depends  on  whether  or  not 
j  =  k,  i.e.,  whether  the  point  is  identical  with  the  point  ^k» 

When  j  k,  then  K^k  is  taken  to  be  the  value  of  the  kernel  function 
K^7,^)  at  ^  and  y  =  £k<  In  particular  the  four  kernel  functions 

of  interest  have  the  following  forms,  with  R  = 


M 


eikR 
jk  “  4*R 


(41) 


ikR 

Djk  =  4rtR~  (ik  '  ^"-k*-^ 


(42) 


ikR  i 

V  -  fcr  <ik  -  V 


(43) 


_  r 

=  4rcR3  L 


(3ik  -  |+  k8R)^R)(Sk-R)  -  (ik  - 


(44) 


When  j  =  k  these  forms  have  singularities  and  are  inapplicable. 
Instead,  in  conformity  with  the  original  definitions  of  these  functions 
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given  in  Appendix  A,  we  assume  that  Kjj  is  a  mean  value  of  the  particular 
kernel  function  over  the  area  surrounding  the  point  ^j>  Furthermore, 
we  assume  that  if  Sj  is  sufficiently  small,  the  mean  value  of  Kjj  will  be 

the  same  as  though  the  surface  element  were  a  spherical  cap  with  £1  at 

2  •* 
its  pole,  with  surface  area  jtb.  =  S.,  and  with  the  same  mean  radius  of 

j  j 

curvature  c^  as  the  actual  surface  at  Then  by  special  calculations 

for  such  a  spherical  cap  (the  procedures  are  similar  to  the  derivations  in 
Appendix  C)  we  obtain  the  following  results  which  ere  correct  to  the  first 
order  in  kb^: 


M  >v  ■  ■ 

JJ  -  2nb . 


ikb,-  1 


DJ]  "  "JJ  *  ^bjCj 


Ejj  ~  2nb3 


Note  that  the  values  of  and  Nj j  are  ordinarily  negligible,  regardless 
of  the  shape  assumed  cor  the  element  of  area,  because  the  curvature  of  the 
element  is  usually  small.  In  fact,  It  is  common  to  approximate  the  neigh* 
borhood  of  each  station  by  a  plane  surface  element  for  which  c^  =  “  and 

djj  '  "jj  ' 

In  summary  then,  the  recommended  procedure  is ;  (1)  pick  N  stations 

on  the  surface;  (ii)  associate  an  area  S.,  a  normal  n.,  and  a  curvature 
•  ^  J  J 

S  with  each  station;  and  (iii)  replace  each  integral  by  a  finite  sum  of 

the  form  of  (40),  where  and  where  is  given  by  (41)  to  (47), 

This  is  probably  the  simplest  and  most  universally  practical  of 
possible  procedures,  but  most  calculators  have  used  more  complex  schemes. 
Usually  they  prescribe  the  shape  of  each  area  element  to  be  a  plane  quadri¬ 
lateral  or  a  triangle  and  then  prescribe  the  location  of  the  station  within 
the  element  of  area.  This  is  equivalent  to  picking  the  weighting  factor 
80  that  it  depends  on  the  shape  of  as  well  as  the  size  of  S^,  The 
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weighting  factor  may  also  depend  on  some  assumed  form  of  variation  of  the 
integrand  in  the  neighborhood  of  These  more  complex  schemes  can  be 

more  accurate  in  special  geometries,  or  they  may  be  more  accurate  when 
calculating  the  main  value  of  each  separate  term  K ^ j •  However,  for  a  body  of 
arbitrary  shape,  it  is  doubtful  that  the  accuracy  of  the  total  calculation 
will  be  improved  at  all.  In  fact,  the  procedure  recommended  here  is  so 
much  more  economical  in  computing  time,  etc.,  that  it  is  often  practical 
to  increase  the  accuracy  by  increasing  the  number  cf  stations. 

One  special  gecmetr\  occurs  so  commonly  in  real  problems,  however, 
as  tc  justify  the  development  of  a  special  calculation  scheme.  These  are 
problems  where  the  vibrating  surface  can  be  approximated  by  a  surface  of 
revolution,  and  where  the  specified  boundary  condition  can  be  written  in 
the  form  of  cos(ncp)  times  an  arbitrary  function  of  the  axial  coordinate, 
where  cp  is  the  angle  about  the  axis  of  revolution.  A  common  and  simple 
example  of  this  problem  would  be  scattering  an  incoming  plane  wave  by  a 
surface  of  revolution.  For  such  a  problem  the  integral  equations  can  all 
be  reduced  to  equations  having  only  one  dimension.  Each  integral  term 
can  be  replaced  by  a  sum  over  a  finite  number  of  stations  in  which  each 
station  is  not  a  point  location  on  the  surface  but  a  circular  ring  on  the 
surface  about  the  axis  of  revolution.  The  elements  of  the  kernel  matrices 
then  have  much  more  complicated  forms  and  are  more  laborious  to  compute, 
although  their  accuracy  is  greatly  increased.  Also  it  is  feasible  to  use 
more  sophisticated  quadrature  formulas,  A  detailed  explanation  of  the 
solution  of  the  Helmholtz  surface  equation  for  this  type  of  problem  is 
given  in  Reference  19. 

Reverting  to  the  general  procedures,  we  find  that  each  of  the 
integral  equations  for  sound  radiation  may  be  replaced  by  a  matrix  equa¬ 
tion.  We  will  denote  a  matrix  by  an  underbar.  Take  £,  £,  C,  and  £  to 
be  column  matrices,  each  with  N  complex  elements  p(^j),  q(£j),  etc;  J>  *s 
a  diagonal  matrix  with  N  diagonal  elements  S j ;  and  M,  D,  N,  and  E  are 
square  matrices  with  N3  complex  elements  defined  in  Equations  (41) 
through  (47),  and  I  is  the  identity  matrix. 

Thus  Equations  (8)  and  (9)  for  the  equivalent  monopole  density 
0(i)  are  replaced  by 
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(2 NS  -  1)0  *  2g 


(48) 


and 

MSP  =  £  (49) 

Equations  (10)  and  (11)  for  the  equivalent  dipole  density  11(2)  are 
replaced  by 


ESu  =  q  (50) 

and 

(2DS  +  I)£  =  2£  (51) 

The  two  independent  integral  equations  (12;  and  (13)  for  the  Helmholtz 
surface  layers  may  be  approximated  by 


(2DS  -  I)£  *  2MSq  (52) 

and 

(2NS  +  1)^  =  2ES£  (53) 

Furthermore,  we  can  form  an  arbitrary  linear  combination  of  (52)  and  (53) 
to  obtain  an  indefinite  number  of  matrix  equations  for  £  or  £. 

Finally  the  Helmholtz  interior  equation  (14)  may  be  replaced  by 


gSq  =  hSp  (54) 

where  £  and  h  are  row  matrices,  not  column  matrices,  with  elements 


gikR 
8x  "  4«R 


(55) 
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with  R  -  x  -  ^  ,  and  x  is  an  arbitrary  fixed  point  within  the  surface. 

If  £  or  £  is  to  be  calculated  completely  from  (54),  there  must  be  at 
least  N  interior  points  x  and  N  equations  such  as  (54). 

A  formal  solution  can  be  written  for  each  of  these  equations  in 
terms  of  the  inverse  matrix.  Thus  the  formal  solutions  to  (48)  and  (49) 
are 


a  =  (2NS  -  I)"1  2£ 

(57) 

£  =  (MS)"1  £ 

'  (58) 

where  the  superscript  -1  denotes  an  inverse  matrix.  However,  when  k  is 
equal  to  a  characteristic  wave  number  of  the  kernel  matrix,  the  matrix 
is  singular,  and  the  inverse  does  not  exist.  The  characteristic  wave 
numbers  of  the  matrix  are  the  same  (except  for  some  shifting  due  to  the 
finite  approximations)  as  the  characteristic  wave  numbers  of  the  kernel 
function  in  the  integral  equation.  The  characteristic  wave  numbers  SL^ 
and  m^  of  the  first  and  second  kind  and  the  characteristic  wave  functions 
pj1^,  q^1',  p^  ,  and  q£“^  which  we  have  previously  defined  in  Equations 
(22)  through  (28)  are  related  to  the  eigenvalues  and  eigenvectors  of  the 
matrices  of  (48)  through  (53).  We  may  consider  that  an  element  Ajnn(k)  of 
a  typical  matrix  A  is  a  function  of  the  wave  number  k  as  well  as  of  the 
two  surface  points  ^  and  and  we  define  the  eigenvalue  'Vj(k)  and  the 
eigenvector  jjj(k)  by  the  equation 


A(k)  «p.(k)  =  v  (k)  j£j(k) 


(59) 


Then,  according  to  (22)  through  (25),  when  k  = 
2MS  has  eigenvalue  0  and  eigenvector  ^ 
2DS  has  eigenvalue  +1  and  eigenvector  jg^1 ^ 
2NS  has  eigenvalue  +1  and  eigenvector 


Also,  according  to  (26)  through  (28),  when  k  = 

(2) 

2DS  has  eigenvalue  -1  and  eigenvector _p^ 

(3) 

2NS  has  eigenvalue  -1  and  eigenvector 

(S) 

2ES  has  eigenvalue  0  and  eigenvector^ 

A  more  complete  analysis  of  the  eigenvalues  and  eigenvectors  of  these 
matrices  for  arbitrary  values  of  the  wave  number  k  will  be  given  in  a 
subsequent  report. 


NUMERICAL  SOLUTIONS  TO  MATRIX  EQUATIONS 


The  matrix  equations  discussed  here  are  all,  except  for  (54),  con¬ 
ventional  linear  matrix  equations  of  the  form  Af  =  £,  There  is  a  vast  and 
growing  literature,  e.g.,  References  (21)  and  (22),  about  numerical 
methods  of  solution,  and  most  of  their  analyses  are  applicable  here.  It 
is  well-known,  for  example,  that  if  is  not  efficient  to  obtain  a  solution 
by  calculating  the  inverse  matrix  A  \  The  methods  of  solutions  are 
commonly  classified  into  direct  methods,  where  the  solution  is  obtained 
after  a  definite  number  of  prescribed  operations,  and  iterative  methods, 
where  the  solution  is  obtained  by  an  indefinite  number  of  successively 
improved  approximations.  In  general  the  choice  of  an  optimum  method  depends 
upon  such  factors  as  the  size,  i.e,,  the  order,  of  the  matrix  A;  whether  the 
matrix  is  sparse,  i.e.,  has  many  zero  elements;  whether  the  matrix  is 
nearly  singular  or  whether  the  problem  is  otherwise  "ill  conditioned";  and 
whether  an  approximate  solution  is  available. 

In  the  particular  problems  of  sound  radiation  and  scattering,  all 
the  matrix  elements  are  finite,  complex,  and  generally  nonzero,  although 
the  elements  of  the  principal  diagonal  of  A  are  generally  dominant.  The 


numbers  of  stations  on  the  surface  often  exceeds  100,  which  means  that 
the  order  of  the  matrix  exceeds  200,  Also,  in  these  problems,  the  wave 


numoer  k  is  often  greater  than  Aj,  the  lowest  characteristic  wave  number 


of  the  first  kind;  in  most  cases  the  characteristic  wave  numbers  and 


m^  in  the  neighborhood  of  k  are  not  known  apriori,  and  it  is  very  laborious 
to  compute  them. 
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When  k  «  ,  then  it  is  a  simple  matter  to  start  with  an  integral 

equation  such  as  (12),  (8),  or  (9)  that  does  not  have  a  singularity  at 
k  =  0.  All  conventional  direct  and  iterative  methods  of  solution  are  then 
applicable,  and  the  optimum  choice  is  mostly  a  matter  of  convenience.  Some¬ 
times  in  this  low-frequency  case,  a  solution  to  the  Laplace  equation,  k  =  0 
in  Equation  (1),  may  be  available  as  an  initial  approximation  in  an  iterative 
solution  to  the  matrix  equation;  however,  when  k  is  not  too  close  to  a 
characteristic  value,  the  problem  is  well  conditioned--even  without  an 
approximate  solution. 

As  already  discussed,  when  the  wave  number  k  is  sufficiently  large, 
k  will  be  arbitrarily  close  to  one  or  more  characteristic  wave  numbers, 
the  matrix  A  will  be  almost  singular,  and  the  problem  will  be  ill-conditioned,, 
Even  when  the  solution  must  simultaneously  satisfy  two  matrix  equations  hav¬ 
ing  no  characteristic  wave  numbers  in  common,  it  appears  likely  that  the 
problem  will  be  ill-conditioned  for  sufficiently  high  k. 

However,  in  the  high-frequency  case,  if  only  the  far-field  is 
ultimately  required,  then  there  is  an  approximation  for  the  surface  pressure 
that  can  be  used  directly  in  (7)  without  solving  any  integral  equation. 

P(l)  (ik  -  -i)  -  q(£)  (60) 

where  c(jr)  is  a  mean  radius  of  curvature  of  the  local  surface  at  £.  Presum¬ 
ably  a  mere  accurate  solution  can  be  obtained  by  using  this  equation  as  the 
initial  approximation  in  an  iterative  solution  of  matrix  equations  (52)  or 
(53)  and  by  terminating  the  iteration  before  the  process  diverges. 

When  the  wave  number  has  some  intermediate  value  greater  than 
but  less  than  the  wave  numbers  at  which  there  is  a  high  density  of  char¬ 
acteristic  values,  then  either  direct  or  iterative  methods  of  solution  may 
be  used.  The  principal  advantage  of  the  direct  methods  is  that  many  library 
programs  are  available  which  can  be  used  automatically.  These  direct 
methods  generally  require  more  storage  space  and  have  larger  computing  time 
(proportional  to  N3)  than  the  iterative  methods  (proportional  to  N3).  The 
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more  common  library  routines  are,  however,  restricted  to  approximately 
100  elements.  A  more  serious  disadvantage  of  a  direct  method  is  that  the 
apparent  solution  could  be  grossly  inaccurate  when  the  wave  number  is  toe 
close  to  a  characteristic  value  without  an  indication  that  the  solution 
is  incorrect.  For  example,  if  an  indeterminate  solution  exists  at  the 
characteristic  value,  then  the  direct  method  will  produce  some  unknown 
multiple  of  the  indeterminate  component;  if  no  finite  solution  exists  at 
the  characteristic  value,  then  the  direct  method  will  tend  to  produce  a 

finite  apparent  solution,  having  an  incorrect  component. 

o' 

A  iterative  method  can  easily  be  programmed  for  any  one  of  these 
matrix  equations.  The  solution  is  ordinarily  faster  than  a  direct  method 
and  requires  a  smaller  storage  space,  so  that  the  number  of  stations 
could  be  increased  with  a  consequent  gain  in  accuracy.  When  the  matrix 
equation  results  from  an  integral  equation  of  the  second  kind,  e„g.  (48), 
it  is  already  in  a  form  to  fit  the  simplest  iteration  scheme,  one  which 
is  the  direct  analog  of  the  classical  Neumann  series  for  solving  the 
integral  equation. 


f(n+1>  =  Bf(n)  +  £  (61) 

where  f^  is  the  n^  approximation  to  f,  and  f^  =  g.  For  example,  in 
(48),  F  =  O,  £  =  -2c[,  and  B  =  2NS.  When  the  matrix  equation  results  from 
an  integral  equation  of  the  first  kind,  e.g.  (49),  it  is  in  the  form 
AF  =  £,  which  can  be  transformed  into  the  form  of  (61)  by  the  scheme 


f(n+1>  =  BF(n>  +  L'1  £ 


(62) 


where  B  =  1^  -  L  1  A,  and  L  is  a  diagonal  matrix  whose  diagonal  elements 
are  equal  to  the  diagonal  elements  of  A.  Note  that  (45)  and  (47)  show 
that  the  diagonal  elements  of  M  or  E  are  real,  finite,  and  of  one  sign 
so  that  the  inverse  matrix  L  1  always  exists  for  matrix  equations  like 
(49),  (50),  and  (52). 
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Equations  (61)  and  (62)  are  called  Jacobi  iteration  schemes,  imply¬ 
ing  that  the  (n+1)  approximation  to  every  element  of  £  is  based  on  the 
(n)*"*1  approximation  to  every  element  of  _f.  In  practice,  some  amount  of 
Causs-Seidel  iteration  is  used,  implying  that  the  (n+1)  approximation  to 
the  j*"*1  element  of  ^  is  based  on  the  (n+1)  approximation  to  the  first 
j-1  elements  of  (which  have  already  been  calculated)  and  to  the  n  approx¬ 

imation  to  the  remaining  elements  of  The  Gauss-Seidel  iteration  usually 
converges  faster  than  does  the  Jacobi  iteration.  The  Gauss-Seidel  scheme 
can  also  be  written  in  the  form  of  (61)  but  with  a  modified  iteration  matrix 
instead  of  B. 

The  principal  disadvantage  of  an  iterative  method  is  tha':  the  itera¬ 
tion  may  not  converge,  even  in  those  cases  where  a  finite  solution  exists. 
Conversely,  an  important  advantage  of  an  iterative  method  is  that  when  a 
finite,  and  possibly  indeterminate,  solution  exists,  and  when  the  iteration 
does  converge,  it  normally  converges  on  the  correct  and  unique  solution 
that  fits  tae  boundary  condition.  An  iterative  process  such  as  (61)  con¬ 
verges  if  and  only  if  the  largest  eigenvalue  of  the  matrix  B,  i.e.,  the 
largest  root  of  the  equation  |b  -  Yl|  =  0,  has  a  modulus  less  than  unity. 
This  condition  is  difficult  to  apply  because  the  eigenvalues  are  difficult 
Co  compute,  unless  the  surface  has  some  special  shape,  e.g.,  a  quadric 
surface.  A  sufficient,  but  not  necessary,  condition  that  the  iteration 
converge  is  that  some  norm  of  B  be  less  than  unity,  ||b| |<  1.  In  our 
limited  experience  with  solutions  of  the  matrix  equations  for  sound  radia¬ 
tion,  a  sufficient  condition  for  the  convergence  of  either  the  Jacobi  or 
Gauss-Seidel  processes  is  that  the  wave  number  k  be  sufficiently  removed 
from  a  characteristic  wave  number  at  which  the  determinant  of  A  =  0.  In 
particular, iterative  solutions  to  (48)  and  (52)  always  seem  to  converge 
as  k  -*  0. 

There  are  also  more  sophisticated  schemes,  which  expand  or  other¬ 
wise  modify  the  range  of  k,  for  which  the  iteration  will  converge  on  a 
correct  and  unique  solution.  In  particular,  there  is  the  iteration  scheme 

f(n+l)  =  hf(n)  +  (1  -  h)  (Bf(n)  +  £)  (63) 


34 


where  h  is  a  complex  number.  In  Reference  17  we  prove  that  this  process 
converges  to  the  correct  solution  provided  all  the  eigenvalues  Yi  of  the 
matrix  B  are  within  a  circle  in  the  complex  plane  passing  through  the 
point  (1,0),  having  a  center  at  h/(h-l),  and  having  a  radius  equal  to  the 
modulus  of  l/(l-h).  When  h  =  0,  this  reduces  to  the  simple  scheme  of  (61). 
When  h  is  a  real  and  positive  number  less  than  unity,  the  domain  of  Y^ 
within  which  the  process  converges  is  a  large  circle  that  includes  the  unit 
circle  of  the  simple  iteration  scheme.  In  fact,  (63)  with  h  =  1/2,  or 
perhaps  a  Gauss-Seidel  variation  of  (63)  with  h  =  1/2,  i c  a  good  wide- 
range  iteration  scheme  that  is  worth  trying  when  the  eigenvalues  Y^  are 
not  known  apriori. 

SPECIAL  METHODS  NEAR  CHARACTERISTIC  WAVE  NUMBERS 

There  probably  is  no  general  method  for  obtaining  the  correct  unique 
solution  when  k  is  close  to  a  k^  for  which  the  solution  is  indeterminate. 

But  there  are  always  special  methods  in  special  situations.  We  have 
already  mentioned  several  such  methods: 

(i)  Shift  the  nearby  k^  by  changing  to  an  integral  equation  with  a 
different  spectrum  of  characteristic  wave  numbers. 

(ii)  Use  the  iteration  scheme  of  Equation  (63)  and  pick  h  30  that  the 
process  converges. 

(iii)  Seek  a  simultaneous  solution  of  two  integral  equations. 

(iv)  Use  the  high-frequency  approximation  of  (60)  and  avoid  the  solution 
of  any  integral,  equation. 

There  are  several  other  methods: 

(v)  The  correct  end  unique  solution  is  clearly  continuous  in  k.  Hence, 
if  the  characteristic  k^  are  not  too  close  together,  we  may  solve  the 
matrix  equation  and  obtain  unique  solutions  with  two  different  values  of 
the  wave  number  on  either  side  of  the  nearest  k^,  and  then  we  can  inter¬ 
polate  to  get  the  correct  solution  at  k^. 

(vi)  Another  method,  which  apparently  has  been  tried  in  electromagnet \c 

23 

wave  theory,  is  to  eliminate  the  possibility  of  internal  standing  waves 
by  adding  an  energy  dissipation  term  to  the  integral  equation.  Presumably, 
if  the  damping  or  dissipation  is  small  enough,  the  solution  to  the  equation 
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with  damping  will  not  differ  appreciably  from  the  correct  solution  without 
damping.  A  simple  way  to  introduce  damping  is  to  treat  k  as  e  complex 
number  with  a  small  imaginary  part.  However,  when  the  real  part  of  k  comes 
too  close  to  a  characteristic  k^,  then  we  can  expect  that  the  matrix  problem 
will  be  ill-conditioned. 

(vii)  There  are  many  situations  where  a  low-frequency  approximation  to 
the  unknown  function  p(y)  tan  be  used  with  negligible  error  in  Equation  (7) 
for  p(^).  The  approximation  might  be  obtained  by  solution  of  the  integral 
equation  with  a  value  of  k  «  Jl^.  Or  the  approximation  might  be  estimated 
from  some  incompressible  potential  flow  (k  a  0)  result  which  does  not 
entail  Ihe  solution  of  any  integral  equation. 

(viii)  As  a  final  example  of  a  special  method  we  solve  the  matrix  equa¬ 
tion  (52)  for  £  by  the  iteration  scheme  of  (63),  and  use  matrix  equation  (53) 
to  determine  a  non- stationary  value  for  the  iteration  parameter  h.  If  £^ 
is  the  n^  approximation  to  £,  and  is  a  tentative  value  for  £^n+1  ^ , 

£(+)  =  2DS£(n)  -  2MS£  (64) 

£(n+l)  -  £(+)  +  h(£(n>  -  £(+))  (65) 

The  parameter  h  is  picked  so  that  £^n+1^  approximately  satisfies  (53)  at 
any  m,  1  £  m  £  N,  stations.  We  truncate  the  N,  1,  and  E  matrices  to 
m  x  N  rectangular  matrices  and  denote  them  by  N^,  and  respectively. 
Then 


S  +  US 


2E  S£^n+1' 


2E 


S£(+)  +h2E 


I  (£ 


(n) 


£(+)) 


(66) 


This  equation  gives  m  different  values  for  the  complex  number  h.  To 

define  a  "best"  value,  we  multiply  on  the  left  by  fS  where  S  is  a 

~ nn  “m 

diagonal  matrix  with  m  diagonal  elements  Sj,  and  f  is  a  row  vector 
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whose  elements  are  the  complex  conjugate-transpose  (hermitian  transposed 

H 

of  the  column  vector  in  the  last  term  of  (66).  That  is,  if  J  denotes 
a  hermitian  transpose 

1H  =  2E  S(£(n)  -  £(+))  (67) 


and 


i  s 


■m 


(2N  S  +  I  )q  =  f  S  2E  Sp(+)  +  h  f  S  fH 
-m  —  — m  - m  — m  — - 


(63) 


and 


f  S  t(2N  S  +  I  )q  -  2E  Sp(+)) 
- m  — m  “  -m  J  nm  — ^ 


f  S  f 
- m  — 


,H 


(69  i 


The  iteration  factor  h  could  be  computed  in  this  way  for  each  itera¬ 
tion,  but  it  seems  safer  to  do  the  first  few  iterations  with  h  =  1/2.  Ihe 
m  stations  for  the  calculation  should  be  selected  at  positions  where  the 
component  of  f,  as  defined  in  (67),  has  a  maximum  modulus.  If  the  itera¬ 
tion  process  has  already  converged  except  for  the  contribution  of  a  single 
characteristic  mode  of  indeterminate  amplitude,  then  in  principle  it  would 
be  sufficient  to  select  only  a  single  station  (m  =  1)  for  the  calculation 
of  h,  but  in  practice  it  would  be  safer  to  use  more  stations. 

In  almost  all  of  these  special  methods,  it  would  be  helpful,  but 
not  necessary,  to  know  the  values  for  the  nearby  characteristic  wave 
numbers  k^,  and  to  have  some  knowledge  of  their  associated  wave  functions 
Pt(j)  and  qi(jy).  For  example,  in  this  last  method,  it  would  be  preferable 
tc  pick  the  stations  for  the  calculation  of  the  iteration  factor  h  at 
these  points  at  which  p^(^)  is  a  maximum. 

SUMMARY  AND  CONCLUSIONS 

We  have  derived  and  discussed  many  different  methods  for  solving 
sound -radiation  and  scattering  problems  associated  with  arbitrary  surfaces. 
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The  common  features  of  all  of  these  methods  are:  (1)  they  replace  the 
radiating  or  scattering  surface  by  an  "equivalent"  surface  distribution 
of  monopoles  and/or  dipoles,  (2)  they  specify  the  density  of  this  equiva¬ 
lent  layer  by  an  integral  equation,  (3)  they  approximate  the  integral 
equation  with  a  matrix  equation  of  finite  order,  (4)  they  solve  the  matrix 
equation  by  the  general  numerical  techniques  which  have  been  developed  for 
matrix  operations  on  high-speed  computers.  We  have  also  derived  and  dis¬ 
cussed  the  conditions  for  the  nonexistence  or  indeterminacy  of  solutions 
to  these  equations,  and  we  have  suggested  several  techniques  for  obtaining 
a  correct  and  unique  solution  when  such  conditions  obtain. 

Some  general  recommendations  have  been  made  for  these  methods  which 
are  likely  to  be  faster,  more  accurate,  or  more  elegant  or  to  have  a  wider 
range  of  application.  However,  it  is  recognized  that  every  one  of  these 
methods  is  complicated  and  laborious  to  implement  and  that  each  method 
may  have  special  advantages  for  special  classes  of  problems.  Furthermore, 
it  happens  that  in  many  real  engineering  problems,  the  physical  parameters 
which  affect  the  solution--e. g. ,  the  details  of  the  boundary  conditions 
or  the  detailed  specification  of  the  surface  shape — can  be  known  only 
approximately  or  only  in  a  statistical  sense.  In  such  cases  there  is 
little  point  in  obtaining  a  .'jolution  whose  mathematical  precision  exceeds 
the  possible  error  due  to  uncertainties  in  the  parameters.  Hence  an 
overall  recommendation  is  to  use  the  simplest  successful  method  that  is 
easily  available  and  tc  give  only  minor  consideration  to  such  factors  as 
speed,  accuracy,  or  wide  range. 

This  has  not  been  a  comprehensive  or  an  historical  review.  For 
example,  we  have  not  covered  what  seems  to  have  been  a  parallel  historical 
development  of  the  same  methods--but  applied  to  problems  in  electromagnetic 
radiation.  A  small  bibliography  on  this  application  is  given  in  Refer¬ 
ence  24.  Also  other  possible  "equivalent"  distribtuions  of  multipoles 

have  not  been  covered,  besides  the  surface  distributions.  For  example,  it 

25 

is  possible  to  find  a  finite  linear  distribution  of  tnonopoles  and/or 

26 

dipoles  inside  a  body  which  solves  arbitrary  Neumann  or  Dirichlet 
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problems  on  the  bodv.  And  it  is  possible  to  express  the  far  field  of  any 

27 

vibrating  body  by  a  suitable  combination  of  point  multipoles.  In  each 
of  these,  equivalent  distributions,  the  multipole  strength  can  probably 
be  specified  by  an  integral  equation  similar  to  those  used  here. 


BLANK 
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APPENDIX  A 


DERIVATION  OF  INTEGRAL  EQUATIONS 

A  common  expression  that  recurs  in  all  the  equations  o  be  derived 
has  the  form  m. /f (R) ,  where  f(R)  is  an  arbitrary  function  of  the  distance 
R,  R  s  £ '  -  z,  m  is  an  arbitrary  unit  vector,  and  Vz  /  means  that  the  gradient 
is  to  be  taken  with  point  z'  varied  and  point  z  fixed.  Then 

m.Vz/f(R)  =  (m.R)  =  -  m.  V,,f(R)  (Al) 


A  y 

where  R  is  a  unit  vector  in  direction  R.  When  the  two  unit  vectors  n 
and  n  are  associated  with  points  z'  and  £.  respectively,  it  is  sometimes 
convenient  to  use  the  abbreviated  notation 


df 


n.V  f 

—  z 


/  df 
dR 


(A2': 


hf 


n'*Vz/f  =  (n<R)  uf 


dR 


(A3; 


To  derive  the  integral  equations  for  the  monopole  layer  O(^),  we 
consider  fir3t  the  field  at  a  point  z',  which  is  a  small  distance  h  off 
surface,  point  in  the  direction  of  the  normal  n/,  i.e.  ,  z  '  -  y_'  t  hn  \ 
We  write  the  field  in  two  parts,  an  integral  over  the  surface  *  within  a 
sphere  of  infinitesimal  radius  6  centered  at  v',  and  an  integral  over 
the  remaining  surface  area  S-A. 


P( 
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ik  z  '  -  ^  | 


4rt  | 


dS  + 


JJ 

S-A 


ae 
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/  1 

£  ‘  1\ 


4jtj  z  '  -  2! 


dS 


;a4) 


Take  R  =  £  -  dS  =  R*dCV(n*R),  where  dft  is  an  element  of  a  solid  angle 

and  evaluate  the  first  integral  as  c  -»  0,  and  subsequently,  as  h  -►  0. 
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dS 


-n 


„  ikR  _s.n  ■.»  _  ikh 

£-10  —ff  21 —  h<m  _  o 

4jtR  (n*R)  C-*oJJ  4«  h-*o 


(A5) 


The  limit  of  the  second  integral  as  e  and  h  approach  zero  in  either  order 
will  be  written  in  an  abbreviated  form  as 


*rr  JS .  ^  jj  £^4-4  as 


J 


**  4n  |  x 7  '  1 


8-40  g_A  Mi'  - 1\ 


(A6) 


Hence,  an  integral  equation  for  a(y) ,  when  p(^)  is  specified  apriori,  is 
given  by 


P(l')  =  *jj  °(l) 


eM  1  'll 

4 Al  m  l\ 


dS 


(A?) 


which  is  the  same  as  (9). 

When  q(^)  is  specified  apriori,  we  proceed  from  (A4)  by  taking  the 
gradient  in  the  n;  direction.  Then 


<1(2  ) 


Aim 
=  /  / 
z  -»y 


r  0  a 

"eM*'  " 

_  4tt!z'  -  2\. 

dS 


f.  If 

z  -*y  JJ 
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The  first  term  becomes 
ikR 


« If ^  <*«■-«  * 


JJ 

A 


_  0(2) 

h-o  2 


(A9) 
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This  term  changes  sign  when  y.1  is  approached  from  the  region  x.  The 
second  integral  in  (A8)  has  the  same  unique  limiting  value,  regardless 
of  the  order  in  which  €  and  h  approach  zero  and  regardless  of  whether 
is  approached  from  the  z  or  x  side.  The  limiting  value  of  this  second 
integr._j.  will  be  written  in  an  abbreviated  notation  as  (with  ft  =  -  2;!) 


(ik  -  p  (n'.R)  dS 


(A10) 


Hence  the  integral  equation  becomes 


which  is  the  same  as  (8). 

To  derive  the  integral  equations  for  the  dipole  density  n(^),  we 
start  with  the  field  at  the  point  z'  =  +  hn'  and  again  separate  the 

form  into  two  parts,  an  integral  over  the  infinitesimal  area  A  about 
and  an  integral  over  the  remaining  surfaces  S-A. 


.  rr  a  Ta11^'  -^1,.  .  rr  a\lkli'-2l 
P<^  >  ’  JJ  11  S:  4,|2'  .  y|  jdS  +  J>  SS  4„|£'  .  ,| 


We  now  follow  the  same  procedures  as  in  (A9)  and  the  first  integral 
becomes 


R3dn 


ikh 

7 -  (1  -  ikh)  dO  ■  —  + 

4tc  n-*o 


') 
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The  limiting  form  of  the  second  integral  in  (A12)  will  be  written  m 
the  abbreviated  form  as 
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(A14) 
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Hence,  the  integral  equation  becomes 
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which  is  the  same  as  (11). 

We  return  now  to  the  field  of  a  dipole  layer  at  the  point 
—  ~  2  +  hn 


pu')  -  ii  i»<x>  « 


(A1 6) 


where  R  =  z  -  We  take  the  gradient  in  the  n/  direction  at  z'  and 
then  take  the  limit  as  z'  -*  y',  assuming  that  the  limit  exists, 


.  /.  Aim 
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dS 
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(A17) 


And  by  repeated  use  of  (A2)  and  (A3), 
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In  fact,  the  gradient  of  the  field  of  a  dipole  layer  _is  continu¬ 
ous  as  the  field  point  crosses  the  dipole  layer.  Thus,  the  limit  of  (A17) 

3 

exists  and  is  unique.  A  rigorous  proof  of  this  result  would  be  long  ana 
tedious  to  develop;  thus,  it  will  be  omitted.  However,  the  result  should 
be  plausible  because  the  dipole  layer  may  be  considered  as  the  superposi¬ 
tion  of  two  monopole  layers,  +  O(y)  =  +  d,  of  equal  amplitude  but 

having  an  opposite  sign  and  being  displaced  an  infinitesimal  distance  d 
apart.  The  gradient  of  the  field  due  to  each  monopole  layer  has  a  dis¬ 
continuity  at  the  surface,  as  given  in  (A9);  however,  in  the  combined 
field  the  two  discontinuities  cancel;  therefore,  the  gradient  of  the  field 
of  a  dipole  layer  has  no  discontinuity  at  the  surface. 

Hence,  an  appropriate  integral  equation  for  tt(^)  when  q(jr)  is 
specified  can  be  stated  in  an  abbreviated  notation  by 


id')  -*n  »(3>  _ 


«■!*'  - i 
dn  [ y  -  yj 


dS 


which  is  the  same  as  (10).  More  explicit  and  detailed  forms  for  the  term 
on  the  right  are  given  by  the  right-hand  sides  of  (A171  or  (A18),  where 
it  must  be  emphasized  that  the  integral  term  in  square  brackets  must  be 
evaluated  with  R  “  \z‘  ~  ^1  before  passing  to  the  limit  z  '  =  y' ,  However, 
if  the  surface  of  integration  is  divided  into  two  parts,  a  small  finite 
area  A  about  y_'  and  the  remaining  surface  S-A,  the  second  part  may  be 
evaluated  with  R  =|^/  -  ,  and  no  limit  i3  necessary  in  the  second  part. 

One  derivation  of  the  Helmholtz  jurface  equation  obtained  by  using 
Green's  transformation  theoram  is  given  in  Reference  15.  However,  once 
we  have  established  (A7)  and  (A15)  for  the  fields  of  arbitrary  monop e 
and  dipole  layers,  we  can  derive  the  Helmholtz  surface  equation  quite 
simply.  We  consider  the  combined  field  p(^)  =  Pm^  +  P<j(2.)  >  where 
Pm(^)  is  the  field  at  the  surface  of  a  monopole  layer  with  surface  density 
minus  q(^) 
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and  pd(y)  is  the  field  at  the  surface  of  a  dipole  layer  with  surface 
density  p(\^)  on  the  positive  side. 


<*'> 


p<^  + 


a  relkh' -ill 


(A21) 


Hence ,  the  combined  field  is 


P(l') 
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(A22) 


which  is  the  same  as  (12). 

The  new  integral  relation  between  the  Helmholtz  surface  layers 
p(^)  and  q(^)  as  given  by  (13)  can  be  derived  by  a  slight  variation  of 
this  last  method.  We  consider  the  combined  normal  gradient  at  the  sur¬ 
face,  q(^')  -  q  (%')  +  qd(^'),  where  q^^7)  is  the  normal  gradient  at 
of  the  field  of  a  layer  of  monopoles  with  surface  density  -q(^),  and 
qd(^7)  is  the  normal  gradient  at  of  the  field  of  a  layer  of  dipoles 
with  surface  density  p(^).  Hence  by  (All)  and  (A17) 
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which  is  the  same  as  (13). 


(A25) 
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APPENDIX  B 


CHARACTERISTIC  WAVE  FUNCTIONS 

The  simple  homogeneous  relations,  (22)  through  (29),  that  apply 
to  the  characteristic  wave  functions  can  all  be  derived  from  the  Green 
transformation  theorems  as  applied  to  an  arbitrary  external  wave  function 
p(z)  with  gradient  q(^)  at  surface  point  %  in  the  direction  of  the  out¬ 
ward  normal  r.. 
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dS 
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The  first  three  equations  are  those  we  have  previously  called  the 
Helmholtz  interior,  surface,  and  exterior  equations  (14),  (12),  and  (7), 
respectively.  Equation  (B4)  is  the  new  surface  gradient  equation  derived 
in  (A25).  Note  that  the  integrals  in  (Bl)  and  (B3)  need  not  be  starred 
because  no  abbreviated  notation  is  implied.  However,  the  possibility  of 
writing  the  set  of  equations  in  this  unified  notation  was  the  principal 
reason  for  selecting  the  particular  form  of  abbreviated  notation. 
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We  also  write  similar  equations  for  any  internal  wave  function 
P(x) ,  with  gradient  Q(^)  in  the  direction  of  the  outward  normal  n. 
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=  P(x) 
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(B5) 
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Note  that  the  forms  of  (El)  through  (B'O  differ  in  sign  from  (B5)  through 
(B8)  only  because  in  the  former  set,  the  normal  points  into  the  region 
where  the  wave  function  p(r)  obeys  the  differential  equation  without 
singularities;  whereas,  in  the  latter  set,  the  normal  points  into  the 
region  where  there  must  be  some  singularities  in  the  wave  function  P(r). 

If  there  were  no  singularities  anywhere,  P(r)  would  be  identically  zero. 

The  special  relations  for  the  characteristic  wave  functions  of 
the  first  kind  would  then  result  from  supplementing  (Bl)  through  (B8) 
with  the  defining  characteristics  for  functions  of  the  first  kind. 

<4°  (2)  -  qjx)(z)  (B9) 
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p<1}Cz)  =  o 


Jb: 


For  example,  (22a)  through  (22c),  respectively,  follow  from  substituting 
(B9)  and  (BlO)  into  (B5)  through  (B7),  respectively. 

The  special  relations  for  the  characteristic  wave  functions  of  the 
second  kind  follow  in  like  manner  from  (Bl)  through  (B8),  plus  the  defin¬ 
ing  characteristics  for  functions  of  the  second  kind 


qf^Cz)  =  o 


We  now  show  that  the  homogeneous  equation 
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has  a  solution  only  when  k  *  L^,  a  characteristic  wave  number  of  the 
first  kind,  i.e.,  only  when  an  internal  wave  P(x)  exists  with  F(^)  =  0. 
For  p(^7)  in  (B13)  can  be  interpreted,  according  to  (All),  as  the  field 
at  £  l'  a  special  dipole  layer  n(^) 


P(Z)  =  *U  H(Z)  J; 
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provided  p(y)  is  taken  equal  to  p(^).  Also,  the  field  of  the  same 
dipole  layer  as  the  point  is  approached  from  the  interior  region  x  is 
given  by 
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In  like  manner,  we  show  that  the  homogeneous  equation 
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( B 1 6  ) 


has  a  solution  only  when  k  =  itk,  i.e.,  only  when  an  interior  standing 
wave  P(x)  exists  with  Q(^)  =  (^P/5n)^  =  0.  First,  note  that  if  a  solu¬ 
tion  to  (B16)  exists  for  a  specific  k,  then  according  to  the  Fredholm 
Alternative  a  solution  must  also  exist  to  the  adjoint  homogeneous 
equation 
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and  by  taking  complex  conjugate  of  (B17),  a  solution  must  alsc  exist 
to  the  homogeneous  equation 


*11  *0 


a  re+ikli'-y| 


33' 


-  il 


dS  = 


Hi) 

2 


(B18) 


Now,  according  to  (8),  £(%  )  can  be  interpreted  as  the  field  gradient 
of  a  special  monopole  layer  o(% ')  as  is  approached  from  z. 
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provided  a(^;)  is  taken  equal  to  However,  the  gradient  at  2 '  °f 

the  internal  field  P(x)  due  to  this  same  monopole  layer  is  given  by 
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Finally,  we  show  that  if  Pi(£)  and  pa (£)  are  two  independent 
solutions  of  the  Helmholtz  equation  with  the  same  wave  number  k,  then 


jJ  Pi  (2)  <*3(2)  dS  “  JJ  Ps(l)  Site)  dS 

Consider  the  Green  formula  as  being  applied  to  the  two  functions  px(jr 
and  pa(z)  in  the  infinite  region  outside  the  surface  y  and  inside  a 
sphere  at  infinity. 
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APPENDIX  C 


SOLUTIONS  TO  THE  INTEGRAL  EQUATIONS 
FOR  THE  PULSATING  SPHERE 

To  solve  the  various  integral  equations  for  the  pulsating  sphere 
it  is  convenient  to  firit  calculate  values  for  the  integrals  of  the  four 
kernel  functions  over  the  surface  of  the  sphere.  We  take  point  z'  to  be 
outside  the  sphere  at  a  distance  r  from  the  origin  at  the  renter.  We 
take  R  and  Rq  as  indicated  in  Figure  1, 


Figure  1  -  Geometry  of  Pulsating  Sphere 
and  take  the  element  of  area  to  be  in  the  shape  of  an  annular  zone  with 

dS  -  2itaasin  0  d0  -  2«a  RdR/r  *  2«R  dR  (Cl' 

0  0  t 

We  can  now  calculate  the  four  integrals  M  ,  D  .  N  ,  and  E  . 

n '  a  '  a  7  a 
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The  integral  equations  (8)  through  (14)  for  this  problem  all  reduce  to 
simple  algebraic  forms  in  terms  of  these  four  values.  For  example,  (13) 
becomes 


qQ(ka  -  tan  ka) 


p0(ik  '  “)  (ka  -  tan  ka> 


(C6) 


or 


q  +  2q  I  “  2p  J 
o  o  o  oo 

which  is  the  same  as  (35).  The  remaining  equations  of  the  series,  (30) 
through  (36),  can  be  obtained  in  a  similar  way. 
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